R Notebook

See https://hbctraining.github.io/Intro-to-ChIPseq/lessons/12_functional_analysis.html

BiocManager::install("ChIPseeker")

'getOption("repos")' replaces Bioconductor standard repositories, see '?repositories' for details

replacement repositories:
    CRAN: https://cran.rstudio.com/

Bioconductor version 3.14 (BiocManager 1.30.16), R 4.1.2 (2021-11-01)
Installing package(s) 'ChIPseeker'
trying URL 'https://bioconductor.org/packages/3.14/bioc/bin/macosx/contrib/4.1/ChIPseeker_1.30.3.tgz'
Content type 'application/x-gzip' length 6936042 bytes (6.6 MB)
==================================================
downloaded 6.6 MB

The downloaded binary packages are in
    /var/folders/c8/rk9y5xz57x5_t_g3v746g7lm0000gn/T//RtmpnOjNYw/downloaded_packages

Old packages: 'ape', 'BH', 'BiocParallel', 'clusterProfiler', 'DBI', 'fs', 'glue', 'igraph', 'knitr', 'Matrix',
  'openssl', 'pkgbuild', 'RcppArmadillo', 'shadowtext', 'tinytex', 'usethis', 'xfun'

Update all/some/none? [a/s/n]: 

a

Warning: unable to access index for repository https://bioconductor.org/packages/3.14/data/annotation/bin/macosx/contrib/4.1:
  cannot open URL 'https://bioconductor.org/packages/3.14/data/annotation/bin/macosx/contrib/4.1/PACKAGES'
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/ape_5.6.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/ape_5.6.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/ape_5.6.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘ape’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/BH_1.78.0-0.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/BH_1.78.0-0.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/BH_1.78.0-0.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘BH’ failed
trying URL 'https://bioconductor.org/packages/3.14/bioc/bin/macosx/contrib/4.1/BiocParallel_1.28.3.tgz'
Content type 'application/x-gzip' length 1375804 bytes (1.3 MB)
==================================================
downloaded 1.3 MB

trying URL 'https://bioconductor.org/packages/3.14/bioc/bin/macosx/contrib/4.1/clusterProfiler_4.2.1.tgz'
Content type 'application/x-gzip' length 711666 bytes (694 KB)
==================================================
downloaded 694 KB

trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/DBI_1.1.2.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/DBI_1.1.2.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/DBI_1.1.2.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘DBI’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/fs_1.5.2.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/fs_1.5.2.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/fs_1.5.2.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘fs’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/glue_1.6.0.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/glue_1.6.0.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/glue_1.6.0.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘glue’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/igraph_1.2.10.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/igraph_1.2.10.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/igraph_1.2.10.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘igraph’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/knitr_1.37.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/knitr_1.37.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/knitr_1.37.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘knitr’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/Matrix_1.4-0.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/Matrix_1.4-0.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/Matrix_1.4-0.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘Matrix’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/openssl_1.4.6.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/openssl_1.4.6.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/openssl_1.4.6.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘openssl’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/pkgbuild_1.3.1.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/pkgbuild_1.3.1.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/pkgbuild_1.3.1.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘pkgbuild’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/RcppArmadillo_0.10.7.5.0.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/RcppArmadillo_0.10.7.5.0.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/RcppArmadillo_0.10.7.5.0.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘RcppArmadillo’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/shadowtext_0.1.0.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/shadowtext_0.1.0.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/shadowtext_0.1.0.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘shadowtext’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/tinytex_0.36.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/tinytex_0.36.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/tinytex_0.36.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘tinytex’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/usethis_2.1.5.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/usethis_2.1.5.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/usethis_2.1.5.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘usethis’ failed
trying URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/xfun_0.29.tgz'
Warning in download.file(url, destfile, method, mode = "wb", ...) :
  URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/xfun_0.29.tgz': status was 'Couldn't resolve host name'
Error in download.file(url, destfile, method, mode = "wb", ...) : 
  cannot open URL 'https://cran.rstudio.com/bin/macosx/contrib/4.1/xfun_0.29.tgz'
Warning in download.packages(pkgs, destdir = tmpd, available = available,  :
  download of package ‘xfun’ failed

The downloaded binary packages are in
    /var/folders/c8/rk9y5xz57x5_t_g3v746g7lm0000gn/T//RtmpnOjNYw/downloaded_packages

library(org.Hs.eg.db)

peakAnno <- annotatePeak(peakfile, tssRegion=c(-3000, 3000), TxDb=txdb)

>> loading peak file...              2021-12-23 20:14:39 
>> preparing features information...         2021-12-23 20:14:39 
>> identifying nearest features...       2021-12-23 20:14:40 
>> calculating distance from peak to TSS...  2021-12-23 20:14:41 
>> assigning genomic annotation...       2021-12-23 20:14:41 
>> assigning chromosome lengths          2021-12-23 20:15:08 
>> done...                   2021-12-23 20:15:08 

peakAnnoDf <- data.frame(peakAnno@anno)
peakEntrez <- peakAnnoDf$geneId
annotations_edb <- AnnotationDbi::select(EnsDb.Hsapiens.v75,
                                         keys = peakEntrez,
                                         columns = c("GENENAME"),
                                         keytype = "ENTREZID")
annotations_edb$ENTREZID <- as.character(annotations_edb$ENTREZID)

ego <- enrichGO(gene = peakEntrez, 
                    keyType = "ENTREZID", 
                    OrgDb = org.Hs.eg.db, 
                    ont = "BP", 
                    pAdjustMethod = "BH", 
                    qvalueCutoff = 0.05, 
                    readable = TRUE)

ekegg <- enrichKEGG(gene = peakEntrez,
                 organism = 'hsa',
                 pvalueCutoff = 0.05)

Reading KEGG annotation online:

Reading KEGG annotation online:

ekegg <- enrichKEGG(gene = peakEntrez,
                 organism = 'hsa',
                 pvalueCutoff = 0.05)
dotplot(ekegg)

promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(peakAnno@anno, windows=promoter)

>> preparing start_site regions by ... 2021-12-22 21:54:07
>> preparing tag matrix...  2021-12-22 21:54:07 

plotAvgProf(tagMatrix, xlim=c(-3000, 3000), conf = 0.95, resample = 1000)

>> plotting figure...            2021-12-22 21:55:08 
>> Running bootstrapping for tag matrix...       2021-12-22 21:55:27 

plotPeakProf2(peak = peakAnno@anno, upstream = rel(0.2), downstream = rel(0.2),
              conf = 0.95, by = "gene", type = "body", nbin = 800,
              TxDb = txdb, weightCol = "V5",ignore_strand = F)

>> binning method is used...2021-12-22 21:55:45
>> preparing body regions by gene... 2021-12-22 21:55:45
>> preparing tag matrix by binning...  2021-12-22 21:55:45 
>> preparing matrix with extension from (TSS-20%)~(TTS+20%)... 2021-12-22 21:55:45
>> 52 peaks(1.028074%), having lengths smaller than 800bp, are filtered... 2021-12-22 21:55:57
>> Running bootstrapping for tag matrix...       2021-12-22 21:59:02 

intron_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = rel(0.2),
                           downstream = rel(0.2), 
                           by = "intron",
                           type = "body",
                           nbin = 100,
                           )

>> binning method is used...2021-12-22 22:08:56
>> preparing body regions by intron... 2021-12-22 22:08:56
>> preparing tag matrix by binning...  2021-12-22 22:08:56 
>> preparing matrix with extension from (TSS-20%)~(TTS+20%)... 2021-12-22 22:08:56
>> 51 peaks(0.1213997%), having lengths smaller than 100bp, are filtered... 2021-12-22 22:09:11

intron_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = rel(0.2),
                           downstream = rel(0.2), 
                           by = "intron",
                           type = "body",
                           nbin = 100,
                           )

>> binning method is used...2021-12-22 22:14:45
>> preparing body regions by intron... 2021-12-22 22:14:45
>> preparing tag matrix by binning...  2021-12-22 22:14:45 
>> preparing matrix with extension from (TSS-20%)~(TTS+20%)... 2021-12-22 22:14:46
>> 51 peaks(0.1213997%), having lengths smaller than 100bp, are filtered... 2021-12-22 22:15:03

plotPeakProf(tagMatrix = intron_matrix, conf = 0.95)

>> Running bootstrapping for tag matrix...       2021-12-22 22:20:50 

intron_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = 1000,
                           downstream = 1000, 
                           by = "intron",
                           type = "start_site",
                           nbin = 250,
                           )

>> binning method is used...2021-12-22 22:32:52
>> preparing start_site regions by intron... 2021-12-22 22:32:52
>> preparing tag matrix by binning...  2021-12-22 22:32:52 

intron_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = 1000,
                           downstream = 1000, 
                           by = "intron",
                           type = "start_site",
                           nbin = 250,
                           )

>> binning method is used...2021-12-22 22:32:57
>> preparing start_site regions by intron... 2021-12-22 22:32:57
>> preparing tag matrix by binning...  2021-12-22 22:32:57 

plotPeakProf(tagMatrix = intron_matrix, conf = 0.95)

>> Running bootstrapping for tag matrix...       2021-12-22 22:33:04 

plotPeakProf(tagMatrix = intron_matrix, conf = 0.95)

>> Running bootstrapping for tag matrix...       2021-12-22 22:47:02 

plotPeakProf(tagMatrix = exon_matrix, conf = 0.95)

>> Running bootstrapping for tag matrix...       2021-12-22 22:44:29 

exon_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = 1000,
                           downstream = 1000, 
                           by = "exon",
                           type = "start_site",
                           nbin=250
                           )

>> binning method is used...2021-12-22 23:31:27
>> preparing start_site regions by exon... 2021-12-22 23:31:27
>> preparing tag matrix by binning...  2021-12-22 23:31:27 

exon_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = 1000,
                           downstream = 1000, 
                           by = "exon",
                           type = "start_site",
                           nbin=250
                           )

>> binning method is used...2021-12-22 23:31:34
>> preparing start_site regions by exon... 2021-12-22 23:31:34
>> preparing tag matrix by binning...  2021-12-22 23:31:34 

plotPeakProf(exon_matrix, xlab="Exon start (5'->3')", ylab = "Peak Frequency", conf = 0.95, resample = 1000)

>> Running bootstrapping for tag matrix...       2021-12-22 23:31:48 

exon_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = 1000,
                           downstream = 1000, 
                           by = "exon",
                           type = "end_site",
                           nbin=250
                           )

>> binning method is used...2021-12-22 23:31:49
>> preparing end_site regions by exon... 2021-12-22 23:31:49
>> preparing tag matrix by binning...  2021-12-22 23:31:49 

plotPeakProf(exon_matrix, xlab="Exon end (5'->3')", ylab = "Peak Frequency", conf = 0.95, resample = 1000)

>> Running bootstrapping for tag matrix...       2021-12-22 23:32:08 

utr_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = 1000,
                           downstream = 1000, 
                           by = "3UTR",
                           type = "start_site",
                           nbin=250
                           )

>> binning method is used...2021-12-22 23:34:17
>> preparing start_site regions by 3UTR... 2021-12-22 23:34:17
>> preparing tag matrix by binning...  2021-12-22 23:34:17 

utr_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = 1000,
                           downstream = 1000, 
                           by = "3UTR",
                           type = "start_site",
                           nbin=250
                           )

>> binning method is used...2021-12-22 23:34:27
>> preparing start_site regions by 3UTR... 2021-12-22 23:34:27
>> preparing tag matrix by binning...  2021-12-22 23:34:27 

plotPeakProf(utr_matrix, xlab="3' UTR start (5'->3')", ylab = "Peak Frequency", conf = 0.95, resample = 1000)

>> Running bootstrapping for tag matrix...       2021-12-22 23:34:33 

utr_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = 1000,
                           downstream = 1000, 
                           by = "3UTR",
                           type = "end_site",
                           nbin=250
                           )

plotPeakProf(utr_matrix, xlab="3' UTR end (5'->3')", ylab = "Peak Frequency", conf = 0.95, resample = 1000)

utr_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = 1000,
                           downstream = 1000, 
                           by = "5UTR",
                           type = "start_site",
                           nbin=250
                           )

>> binning method is used...2021-12-22 23:36:22
>> preparing start_site regions by 5UTR... 2021-12-22 23:36:22
>> preparing tag matrix by binning...  2021-12-22 23:36:22 

utr_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = 1000,
                           downstream = 1000, 
                           by = "5UTR",
                           type = "start_site",
                           nbin=250
                           )

>> binning method is used...2021-12-22 23:36:29
>> preparing start_site regions by 5UTR... 2021-12-22 23:36:29
>> preparing tag matrix by binning...  2021-12-22 23:36:29 

plotPeakProf(utr_matrix, xlab="5' UTR start (5'->3')", ylab = "Peak Frequency", conf = 0.95, resample = 1000)

>> Running bootstrapping for tag matrix...       2021-12-22 23:36:35 

utr_matrix <- getTagMatrix(peak = peakAnno@anno, 
                           TxDb = txdb,
                           upstream = 1000,
                           downstream = 1000, 
                           by = "5UTR",
                           type = "end_site",
                           nbin=250
                           )

plotPeakProf(utr_matrix, xlab="5' UTR end (5'->3')", ylab = "Peak Frequency", conf = 0.95, resample = 1000)

plotPeakProf(gene_matrix, xlab="Gene start (5'->3')", ylab = "Peak Frequency", conf = 0.95, resample = 1000)

>> Running bootstrapping for tag matrix...       2021-12-23 20:16:29 

plotPeakProf(gene_matrix, xlab="Gene end (5'->3')", ylab = "Peak Frequency", conf = 0.95, resample = 1000)

>> Running bootstrapping for tag matrix...       2021-12-23 20:16:06